/* 
 * File:   random01.cpp
 * Author: karol
 * 
 * Created on 17 listopad 2009, 12:50
 */

#include "Random.h"
#ifndef OMP_DEBUG
#include <omp.h>
#else
#include "omp_debug.h"
#endif
/** Definition of RNG object **/

//rangen_wd random01;
//rangen random01;
//rangen_mwc random01;
rangen_mwc random01;


/** plusminusone() */
int plusminusone()
{
    if(random01()>0.5)
        return 1;
    else return -1;
}

/** mwc_wrap **/
std::vector<rangen_mwc> mwc_wrap::rg;
void mwc_wrap::setup(int n)
{
    rg.resize(n);
    for(int i=0;i<n;i++)
        rg[i].setup(i);
}
double mwc_wrap::operator()()
{
    //static std::vector<rng> rg(omp_get_num_threads());
    return rg[omp_get_thread_num()]();
}

/** WydroRNG **/
long int seed1 = 46723;			// seed for random long integer generator
long int seed2 = 2593459;		// seed for shuffling algorithm
long int jj =89343;
double denom=(1.0/(m-1));
const int NR=64;
long int y;			// random int mixed into shuffling array.
long int j[NR];		// vector of random integers that are shuffled
					// Both y and j[] are seeded by making the singl
long int random_int(long int& j)
{
	long int l;
	l = j/q;				// Parameter "j" is recursively defined to produce a
							// a sequence of random numbers.
	j = a*(j - q*l) -r*l;	// This iteration step is actually a mod 2^31
							// operation as described by Schrage.  Here, first term =
							// j mod q and second term is (int divide(j/q) )*r.
							// See Barkema, page 388.
	if (j<0) j += m;		// This shift bring the number back into the interval
							// [1, m - 1].
	return (j-1);	// j - 1 is a random integer in interval [0, m - 2]
}	//end program to return random long integer

double drandom(long int& seed2)
{
	long int l;
	long int k;

	l = seed2/q;
	seed2 = a*(seed2 - q*l) - r*l;	//	Selects "seed2" to be a random integer in
							// [1, 2^31 - 2 ] and updates "seed2" as reference variable.
							//	This ensures that a series of calls will produce a
							//	a series of random long ints "seed2".
	if (seed2 < 0)			//ensures that seed2 is random int in [1, m - 1]
		{
			seed2 = seed2 + m;
		}
	k = (y/m)*NR;		//Selects an integer in [0, NR]

	y = j[k];			// Sets y = integer of k-th component of j[].
						// Side effect: global variable y is updated for
						// next call of drandom.
	j[k] = seed2;		// Side effect: global variable j[k] is updated
						// for next call of drandom.

	return denom*(y - 1);	// Returns randomly shuffled double float in [0, 1).
}	//end Shuffler program

void initializeRNG(long int& seed1)
{
	int s;
	for(s = 0; s< NR ; s++)	//loop iteratively seeds each component of j[].
		{
			j[s] = random_int(seed1); //seeds component "s" of shuffler vector j[].
			//	std::cout << "j[" << s << "] = " << j[s] << '\n';
			//	The above line is used as a test to insure that j[] is a random array
			//	That is, to print out array j[] to check random initialization
		}
	y = random_int(seed1);	// seeds variable j
	//	std::cout << "yy = " << yy<< '\n';
	//	That is, to print out yy to check random initialization
}	//end program to initialize Shuffler program


/** rangen_mwc **/
rangen_mwc::rangen_mwc()
{
    setup(omp_get_thread_num());
}

rangen_mwc::rangen_mwc(const rangen_mwc &)
{
    setup(omp_get_thread_num());
}
void rangen_mwc::setup(int k)
{
    lsp[0][0]=489UL;
    lsp[0][1]=2100239007743UL;
    lsp[0][2]=1050119503871UL;
    lsp[1][0]=1170UL;
    lsp[1][1]=5025111736319UL;
    lsp[1][2]=2512555868159UL;
    lsp[2][0]=1245UL;
    lsp[2][1]=5347234283519UL;
    lsp[2][2]=2673617141759UL;
    lsp[3][0]=1440UL;
    lsp[3][1]=6184752906239UL;
    lsp[3][2]=3092376453119UL;
    lsp[4][0]=1728UL;
    lsp[4][1]=7421703487487UL;
    lsp[4][2]=3710851743743UL;
    lsp[5][0]=2094UL;
    lsp[5][1]=8993661517823UL;
    lsp[5][2]=4496830758911UL;
    lsp[6][0]=2100UL;
    lsp[6][1]=9019431321599UL;
    lsp[6][2]=4509715660799UL;
    lsp[7][0]=2355UL;
    lsp[7][1]=10114647982079UL;
    lsp[7][2]=5057323991039UL;
    lsp[8][0]=2598UL;
    lsp[8][1]=11158325035007UL;
    lsp[8][2]=5579162517503UL;
    lsp[9][0]=2793UL;
    lsp[9][1]=11995843657727UL;
    lsp[9][2]=5997921828863UL;
    lsp[10][0]=2988UL;
    lsp[10][1]=12833362280447UL;
    lsp[10][2]=6416681140223UL;
    lsp[11][0]=3195UL;
    lsp[11][1]=13722420510719UL;
    lsp[11][2]=6861210255359UL;
    lsp[12][0]=3213UL;
    lsp[12][1]=13799729922047UL;
    lsp[12][2]=6899864961023UL;
    lsp[13][0]=3849UL;
    lsp[13][1]=16531329122303UL;
    lsp[13][2]=8265664561151UL;
    lsp[14][0]=3939UL;
    lsp[14][1]=16917876178943UL;
    lsp[14][2]=8458938089471UL;
    lsp[15][0]=4005UL;
    lsp[15][1]=17201344020479UL;
    lsp[15][2]=8600672010239UL;
    lsp[16][0]=4173UL;
    lsp[16][1]=17922898526207UL;
    lsp[16][2]=8961449263103UL;
    lsp[17][0]=4299UL;
    lsp[17][1]=18464064405503UL;
    lsp[17][2]=9232032202751UL;
    lsp[18][0]=4815UL;
    lsp[18][1]=20680267530239UL;
    lsp[18][2]=10340133765119UL;
    lsp[19][0]=4833UL;
    lsp[19][1]=20757576941567UL;
    lsp[19][2]=10378788470783UL;
    lsp[20][0]=5439UL;
    lsp[20][1]=23360327122943UL;
    lsp[20][2]=11680163561471UL;
    lsp[21][0]=5538UL;
    lsp[21][1]=23785528885247UL;
    lsp[21][2]=11892764442623UL;
    lsp[22][0]=5859UL;
    lsp[22][1]=25164213387263UL;
    lsp[22][2]=12582106693631UL;
    lsp[23][0]=6849UL;
    lsp[23][1]=29416231010303UL;
    lsp[23][2]=14708115505151UL;
    lsp[24][0]=7713UL;
    lsp[24][1]=33127082754047UL;
    lsp[24][2]=16563541377023UL;
    lsp[25][0]=9075UL;
    lsp[25][1]=38976828211199UL;
    lsp[25][2]=19488414105599UL;
    lsp[26][0]=9150UL;
    lsp[26][1]=39298950758399UL;
    lsp[26][2]=19649475379199UL;
    lsp[27][0]=9264UL;
    lsp[27][1]=39788577030143UL;
    lsp[27][2]=19894288515071UL;
    lsp[28][0]=10158UL;
    lsp[28][1]=43628277792767UL;
    lsp[28][2]=21814138896383UL;
    lsp[29][0]=10209UL;
    lsp[29][1]=43847321124863UL;
    lsp[29][2]=21923660562431UL;
    lsp[30][0]=10395UL;
    lsp[30][1]=44646185041919UL;
    lsp[30][2]=22323092520959UL;
    lsp[31][0]=10410UL;
    lsp[31][1]=44710609551359UL;
    lsp[31][2]=22355304775679UL;
    lsp[32][0]=10638UL;
    lsp[32][1]=45689862094847UL;
    lsp[32][2]=22844931047423UL;
    lsp[33][0]=10869UL;
    lsp[33][1]=46681999540223UL;
    lsp[33][2]=23340999770111UL;
    lsp[34][0]=10974UL;
    lsp[34][1]=47132971106303UL;
    lsp[34][2]=23566485553151UL;
    lsp[35][0]=11499UL;
    lsp[35][1]=49387828936703UL;
    lsp[35][2]=24693914468351UL;
    lsp[36][0]=12489UL;
    lsp[36][1]=53639846559743UL;
    lsp[36][2]=26819923279871UL;
    lsp[37][0]=12984UL;
    lsp[37][1]=55765855371263UL;
    lsp[37][2]=27882927685631UL;
    lsp[38][0]=13200UL;
    lsp[38][1]=56693568307199UL;
    lsp[38][2]=28346784153599UL;
    lsp[39][0]=14148UL;
    lsp[39][1]=60765197303807UL;
    lsp[39][2]=30382598651903UL;
    lsp[40][0]=14430UL;
    lsp[40][1]=61976378081279UL;
    lsp[40][2]=30988189040639UL;
    lsp[41][0]=14895UL;
    lsp[41][1]=63973537873919UL;
    lsp[41][2]=31986768936959UL;
    lsp[42][0]=15189UL;
    lsp[42][1]=65236258258943UL;
    lsp[42][2]=32618129129471UL;
    lsp[43][0]=15405UL;
    lsp[43][1]=66163971194879UL;
    lsp[43][2]=33081985597439UL;
    lsp[44][0]=15594UL;
    lsp[44][1]=66975720013823UL;
    lsp[44][2]=33487860006911UL;
    lsp[45][0]=15870UL;
    lsp[45][1]=68161130987519UL;
    lsp[45][2]=34080565493759UL;
    lsp[46][0]=16170UL;
    lsp[46][1]=69449621176319UL;
    lsp[46][2]=34724810588159UL;
    lsp[47][0]=16728UL;
    lsp[47][1]=71846212927487UL;
    lsp[47][2]=35923106463743UL;
    lsp[48][0]=17193UL;
    lsp[48][1]=73843372720127UL;
    lsp[48][2]=36921686360063UL;
    lsp[49][0]=17370UL;
    lsp[49][1]=74603581931519UL;
    lsp[49][2]=37301790965759UL;
    lsp[50][0]=18114UL;
    lsp[50][1]=77799037599743UL;
    lsp[50][2]=38899518799871UL;
    lsp[51][0]=18459UL;
    lsp[51][1]=79280801316863UL;
    lsp[51][2]=39640400658431UL;
    lsp[52][0]=18705UL;
    lsp[52][1]=80337363271679UL;
    lsp[52][2]=40168681635839UL;
    lsp[53][0]=19164UL;
    lsp[53][1]=82308753260543UL;
    lsp[53][2]=41154376630271UL;
    lsp[54][0]=19173UL;
    lsp[54][1]=82347407966207UL;
    lsp[54][2]=41173703983103UL;
    lsp[55][0]=19530UL;
    lsp[55][1]=83880711290879UL;
    lsp[55][2]=41940355645439UL;
    lsp[56][0]=19593UL;
    lsp[56][1]=84151294230527UL;
    lsp[56][2]=42075647115263UL;
    lsp[57][0]=19725UL;
    lsp[57][1]=84718229913599UL;
    lsp[57][2]=42359114956799UL;
    lsp[58][0]=19788UL;
    lsp[58][1]=84988812853247UL;
    lsp[58][2]=42494406426623UL;
    lsp[59][0]=19953UL;
    lsp[59][1]=85697482457087UL;
    lsp[59][2]=42848741228543UL;
    lsp[60][0]=21108UL;
    lsp[60][1]=90658169683967UL;
    lsp[60][2]=45329084841983UL;
    lsp[61][0]=21549UL;
    lsp[61][1]=92552250261503UL;
    lsp[61][2]=46276125130751UL;
    lsp[62][0]=22020UL;
    lsp[62][1]=94575179857919UL;
    lsp[62][2]=47287589928959UL;
    lsp[63][0]=22065UL;
    lsp[63][1]=94768453386239UL;
    lsp[63][2]=47384226693119UL;
    lsp[64][0]=22095UL;
    lsp[64][1]=94897302405119UL;
    lsp[64][2]=47448651202559UL;
    lsp[65][0]=22098UL;
    lsp[65][1]=94910187307007UL;
    lsp[65][2]=47455093653503UL;
    lsp[66][0]=22209UL;
    lsp[66][1]=95386928676863UL;
    lsp[66][2]=47693464338431UL;
    lsp[67][0]=22419UL;
    lsp[67][1]=96288871809023UL;
    lsp[67][2]=48144435904511UL;
    lsp[68][0]=22533UL;
    lsp[68][1]=96778498080767UL;
    lsp[68][2]=48389249040383UL;
    lsp[69][0]=22923UL;
    lsp[69][1]=98453535326207UL;
    lsp[69][2]=49226767663103UL;
    lsp[70][0]=23184UL;
    lsp[70][1]=99574521790463UL;
    lsp[70][2]=49787260895231UL;
    lsp[71][0]=23355UL;
    lsp[71][1]=100308961198079UL;
    lsp[71][2]=50154480599039UL;
    lsp[72][0]=24024UL;
    lsp[72][1]=103182294319103UL;
    lsp[72][2]=51591147159551UL;
    lsp[73][0]=24078UL;
    lsp[73][1]=103414222553087UL;
    lsp[73][2]=51707111276543UL;
    lsp[74][0]=24513UL;
    lsp[74][1]=105282533326847UL;
    lsp[74][2]=52641266663423UL;
    lsp[75][0]=24708UL;
    lsp[75][1]=106120051949567UL;
    lsp[75][2]=53060025974783UL;
    lsp[76][0]=25404UL;
    lsp[76][1]=109109349187583UL;
    lsp[76][2]=54554674593791UL;
    lsp[77][0]=25494UL;
    lsp[77][1]=109495896244223UL;
    lsp[77][2]=54747948122111UL;
    lsp[78][0]=26298UL;
    lsp[78][1]=112949049950207UL;
    lsp[78][2]=56474524975103UL;
    lsp[79][0]=26400UL;
    lsp[79][1]=113387136614399UL;
    lsp[79][2]=56693568307199UL;
    lsp[80][0]=26790UL;
    lsp[80][1]=115062173859839UL;
    lsp[80][2]=57531086929919UL;
    lsp[81][0]=28560UL;
    lsp[81][1]=122664265973759UL;
    lsp[81][2]=61332132986879UL;
    lsp[82][0]=28860UL;
    lsp[82][1]=123952756162559UL;
    lsp[82][2]=61976378081279UL;
    lsp[83][0]=29169UL;
    lsp[83][1]=125279901057023UL;
    lsp[83][2]=62639950528511UL;
    lsp[84][0]=29424UL;
    lsp[84][1]=126375117717503UL;
    lsp[84][2]=63187558858751UL;
    lsp[85][0]=29454UL;
    lsp[85][1]=126503966736383UL;
    lsp[85][2]=63251983368191UL;
    lsp[86][0]=29709UL;
    lsp[86][1]=127599183396863UL;
    lsp[86][2]=63799591698431UL;
    lsp[87][0]=30843UL;
    lsp[87][1]=132469676310527UL;
    lsp[87][2]=66234838155263UL;
    lsp[88][0]=30903UL;
    lsp[88][1]=132727374348287UL;
    lsp[88][2]=66363687174143UL;
    lsp[89][0]=31794UL;
    lsp[89][1]=136554190209023UL;
    lsp[89][2]=68277095104511UL;
    lsp[90][0]=31989UL;
    lsp[90][1]=137391708831743UL;
    lsp[90][2]=68695854415871UL;
    lsp[91][0]=32208UL;
    lsp[91][1]=138332306669567UL;
    lsp[91][2]=69166153334783UL;
    lsp[92][0]=32235UL;
    lsp[92][1]=138448270786559UL;
    lsp[92][2]=69224135393279UL;
    lsp[93][0]=32634UL;
    lsp[93][1]=140161962737663UL;
    lsp[93][2]=70080981368831UL;
    lsp[94][0]=32700UL;
    lsp[94][1]=140445430579199UL;
    lsp[94][2]=70222715289599UL;
    lsp[95][0]=33213UL;
    lsp[95][1]=142648748802047UL;
    lsp[95][2]=71324374401023UL;
    lsp[96][0]=33540UL;
    lsp[96][1]=144053203107839UL;
    lsp[96][2]=72026601553919UL;
    lsp[97][0]=34758UL;
    lsp[97][1]=149284473274367UL;
    lsp[97][2]=74642236637183UL;
    lsp[98][0]=35160UL;
    lsp[98][1]=151011050127359UL;
    lsp[98][2]=75505525063679UL;
    lsp[99][0]=35199UL;
    lsp[99][1]=151178553851903UL;
    lsp[99][2]=75589276925951UL;
    lsp[100][0]=35370UL;
    lsp[100][1]=151912993259519UL;
    lsp[100][2]=75956496629759UL;
    lsp[101][0]=36384UL;
    lsp[101][1]=156268090097663UL;
    lsp[101][2]=78134045048831UL;
    lsp[102][0]=36693UL;
    lsp[102][1]=157595234992127UL;
    lsp[102][2]=78797617496063UL;
    lsp[103][0]=36888UL;
    lsp[103][1]=158432753614847UL;
    lsp[103][2]=79216376807423UL;
    lsp[104][0]=37005UL;
    lsp[104][1]=158935264788479UL;
    lsp[104][2]=79467632394239UL;
    lsp[105][0]=37110UL;
    lsp[105][1]=159386236354559UL;
    lsp[105][2]=79693118177279UL;
    lsp[106][0]=37344UL;
    lsp[106][1]=160391258701823UL;
    lsp[106][2]=80195629350911UL;
    lsp[107][0]=38760UL;
    lsp[107][1]=166472932392959UL;
    lsp[107][2]=83236466196479UL;
    lsp[108][0]=38835UL;
    lsp[108][1]=166795054940159UL;
    lsp[108][2]=83397527470079UL;
    lsp[109][0]=39420UL;
    lsp[109][1]=169307610808319UL;
    lsp[109][2]=84653805404159UL;
    lsp[110][0]=39444UL;
    lsp[110][1]=169410690023423UL;
    lsp[110][2]=84705345011711UL;
    lsp[111][0]=40005UL;
    lsp[111][1]=171820166676479UL;
    lsp[111][2]=85910083338239UL;
    lsp[112][0]=40113UL;
    lsp[112][1]=172284023144447UL;
    lsp[112][2]=86142011572223UL;
    lsp[113][0]=40284UL;
    lsp[113][1]=173018462552063UL;
    lsp[113][2]=86509231276031UL;
    lsp[114][0]=40500UL;
    lsp[114][1]=173946175487999UL;
    lsp[114][2]=86973087743999UL;
    lsp[115][0]=40845UL;
    lsp[115][1]=175427939205119UL;
    lsp[115][2]=87713969602559UL;
    lsp[116][0]=41313UL;
    lsp[116][1]=177437983899647UL;
    lsp[116][2]=88718991949823UL;
    lsp[117][0]=41574UL;
    lsp[117][1]=178558970363903UL;
    lsp[117][2]=89279485181951UL;
    lsp[118][0]=41973UL;
    lsp[118][1]=180272662315007UL;
    lsp[118][2]=90136331157503UL;
    lsp[119][0]=42285UL;
    lsp[119][1]=181612692111359UL;
    lsp[119][2]=90806346055679UL;
    lsp[120][0]=42588UL;
    lsp[120][1]=182914067202047UL;
    lsp[120][2]=91457033601023UL;
    lsp[121][0]=42735UL;
    lsp[121][1]=183545427394559UL;
    lsp[121][2]=91772713697279UL;
    lsp[122][0]=43050UL;
    lsp[122][1]=184898342092799UL;
    lsp[122][2]=92449171046399UL;
    lsp[123][0]=43119UL;
    lsp[123][1]=185194694836223UL;
    lsp[123][2]=92597347418111UL;
    lsp[124][0]=43638UL;
    lsp[124][1]=187423782862847UL;
    lsp[124][2]=93711891431423UL;
    lsp[125][0]=45585UL;
    lsp[125][1]=195786084188159UL;
    lsp[125][2]=97893042094079UL;
    lsp[126][0]=45720UL;
    lsp[126][1]=196365904773119UL;
    lsp[126][2]=98182952386559UL;
    lsp[127][0]=45729UL;
    lsp[127][1]=196404559478783UL;
    lsp[127][2]=98202279739391UL;

    unsigned long seed = lsp[k][2];
    unsigned long mult = lsp[k][0];

    gen.RandomInit(seed,mult);
}
double rangen_mwc::operator()()
{
    return gen.Random();
}

/** rangen **/
boost::rand48 rng2;
rangen::rangen():rng(rng2()),uni01(rng){}
rangen::rangen(const rangen &):rng(rng2()),uni01(rng)
{
}
double rangen::Gen01()
{
    return uni01();
}
double rangen::operator()()
{
    return uni01();
}

/** rangen_wd **/
rangen_wd::rangen_wd()
{
    initializeRNG(seed1);
}
double rangen_wd::Gen01()
{
    return drandom(seed2);
}
double rangen_wd::operator()()
{
    return Gen01();
}
